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We report on probability-density-functions (PDF) of the mass density in numerical simulations of 
highly compressible hydrodynamic flows and the corresponding structure formation of Lagrangian 
particles advected by the flows. Numerical simulations were performed with 512'^ collocation points 
and 2 million tracer particles integrated over several dynamical times. We propose a connection 
between the PDF of the Lagrangian tracer particles and the predicted log-normal distribution of 
the density fluctuations in isothermal systems. 
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Introduction: Compressible fluid and plasma turbulence 
is believed to be of fundamental importance for under- 
standing relevant processes in astrophysics. A typical 
example is the effect of turbulence on the star formation 
rate in dense molecular clouds in the interstellar medium 
[H . The effect of turbulence is manifested in the highly 
intermittent nature of density fluctuations in supersonic 
flows which naturally has an enormous impact on the 
local conditions triggering star formation. 

Traditionally, turbulence is studied in the Eulcrian de- 
scription, where the temporal evolution of the fluid quan- 
tities is followed at flxed spatial locations. However, 
in the last 10 years the Lagrangian description, which 
describes the motion of fluid particles, has undergone 
a rapid development. This is mainly due to advanced 
experimental techniques introduced in Ris0 0, Cornell 
0, 0, 0] and Lyon p for studying the statistics of pas- 
sive tracer particles in incompressible turbulent fluids. 
Both, experimental and numerical 0, Q studies highlight 
the importance of singular structures like vortex tubes 
and sheets for the intermittent statistics of incompress- 
ible flows. Lagrangian statistics is not only interesting 
for obtaining a deeper understanding of the influence of 
typical coherent or nearly-singular structures in the flow 
but also of fundamental importance for understanding 
mixing, clustering and diffusion properties of turbulent 
astrophysical fluid and plasma flows. 

In this Letter, we will use the Lagrangian description 
of passive tracer particles in highly compressible flows. 
Due to the compressibility, the situation is very different 
from the incompressible case. Clustering of passive tracer 
particles will reflect the strong density fluctuations. This 
effect does not appear in the incompressible case. How- 
ever, the incompressible case will serve as a guideline to 
obtain a mapping between the PDF of density fluctua- 
tions and the PDF of the particle distribution. 

The goal of the paper is threefold: The main result is 
an analytical expression for the spatial distribution of the 
Lagrangian tracers. This result will be used to confirm 
and quantify the relation between the standard deviation 
of the Eulerian density fluctuations and the mean Mach 
number. And last, the analytical description helps to 



define a precise structuring time for the clustering of the 
Lagrangian particles which turns out to be significantly 
shorter than the dynamical time. 

Numerical Methods: To follow this approach in the ex- 
amination of compressible fluid turbulence we solved the 
isothermal Euler equations for density p and momentum 
density u = pv 

dtp + S/-u = 

9tu + V- ( — ) = -VP + pk (1) 
P 

with P ^ p using direct numerical simulations in real 
space, and realized Lagrangian measurements with tracer 
particles, which are advected by the simulated flow. We 
performed three runs of stationary flows (a, b and c) 
with mean r.m.s. Mach numbers Af of 0.4, 1.4 and 4.6, 
respectively. 

Equations ([T]) were solved in the racoon environment 
for hyperbolic differential equations 0. It utilizes a 
third order semi-discrete central scheme [lo| . Integra- 
tion is done by a third order strongly stable Runge-Kutta 
scheme . A 3-dimensional cube of length L = 27: was 
discretized in physical space with = 512^ collocation 
points and periodic boundary conditions. The simula- 
tion was started with an Orszag-Tang like initial condi- 
tion and a homogeneous density. The simulations were 
run for several dynamical times tdyn = L/2M to let the 
system evolve into a turbulent state. The mean energy 
input was realized by a temporal delta correlated stirring 
force with zero mean momentum. Forcing was applied in 
a shell /c e [1,2] in Fourier space, such that it acts only 
on scales comparable to the entire box. The amplitudes 
of the Fourier modes were adjusted in order to achieve 
stationarity at various mean Mach numbers. Then the 
mean energy input was kept constant for each run. 

In this evolved turbulent system we inserted 
N = 2 000 000 particles. They were randomly distributed 
over the computational domain and then advected by 
the flow field v(x) = Particles were integrated us- 

ing the same third-order Runga-Kutta scheme as for the 
fields. The field values at the particle positions x were 
evaluated by a linear interpolation, since the density has 
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FIG. 1: Compensated velocity spectra 



to be positive even in regions where strong shocks cause 
steep gradients close to flat areas. Higher order interpo- 
lators here tend to undershoot which leads to artificial 
particle reflection. After another dynamical time tdyn to 
assure that the particle statistics has reached a station- 
ary state, data were analyzed both for Eulcrian and the 
Lagrangian statistics. 

Eulerian statistics: To characterize the flows we mea- 
sured some specific values and show velocity-spectra and 
mass-density PDFs in the turbulent regime (see Figs. [U 
[2] and Table |T|. In fully developed isotropic turbulence, 
the energy brought into the system at large scales should 
cascade through the scales and finally dissipate at small 
scales, known as the Richardson cascade of energy. In 
stationary turbulent fiows, the energy injection equals 
the dissipation. Therefore, the mean energy dissipation e 
can be obtained from the mean energy input (dtEk) of the 
driving force k averaged over several dynamical times. 
Normalized to the mean kinetic energy in the box £/-Ekin 
we see in Table U that dissipation rises with the Mach 
number. In order to quantify the Eulerian compressibil- 
ity the dimensionless ratio C = ((9^?;^)^) / is 
measured which is zero for solenoidal and one for po- 
tential fiows. This quantity also shows a dependence of 
Mach number. Finally we give an estimate for the inte- 
gral scale Reynolds number R = IqVq/u defined by the 
integral length scale Iq — /e. Here, vq is the 

root mean square velocity fluctuation, which in isother- 
mal flows is the same as the mean Mach number, and 
the dynamical viscosity v is given hy v — jfl'^e^l'^. The 
Kolmogorov dissipation scale 77 was estimated by the dis- 
cretization scale T] = Ax = L/Af. Please note, that this 
is a conservative estimate since in spectral simulations of 
incompressible flows the dissipation scale rj is normally 
half the grid spacing. Only for high Mach number flows 
a Reynolds number can be reached which is comparable 
to the one obtained from spectral simulations of incom- 
pressible flows with the same resolution. The power spec- 



FIG. 2; PDF of logarithm of mass density. Solid lines are 
normal distribution with mean and variance obtained from 
data. 



trum of the velocity fluctuations steepens with increas- 
ing Reynolds number reflecting the occurrence of strong 
shocks. The slope of the spectrum in the inertial range 
drops to a value, which is compatible to P{k) ~ fe"-'^-*^ 
for run c. This is in agreement with the simulations of 
Kritsuk et al. [12] and Boldyrev et al. [l^. Although the 
latter is obtained from MHD simulations, it just reveals 
the dominance of strong shocks that converge to Burgers- 
like spectra ~ for even higher Mach-numbers. 

In the isothermal case the mass density PDF R{p) 
should follow a log-normal distribution as discussed in 
[H . [lit . In Fig. [2] we compared the density PDF of 
the numerical simulations with normal distributions with 
mean and variance obtained from the numerical data 



which confirms the reasoning of IJ, [15 1 . 
Distribution of tracer particles: To obtain the PDF of 
the spatial distribution of the tracer particles we divided 
the computational box into n = 40'^ equally spaced sub- 
domains and counted the particles in them. The number 
of boxes with k particles gave the PDF of the tracers 
as a discrete distribution T{k). Since the particles were 
initially randomly distributed, their PDF started from a 
Poissonian 



Tinit {k) 



with A = N/n and changed in time towards a stationary 
distribution. 



TABLE I: Characteristic quantities for the three runs 



Run 


Mach 


C 




Reynolds number 


a 


0.4 


0.07 


0.31 


250 


b 


1.4 


0.26 


0.47 


700 


c 


4.6 


0.50 


0.95 


1756 



3 



To understand the relation between the mass density 
PDF and the particle distribution, we first consider an 
incompressible flow. The PDF of p (x) ~ po = const is 



Ri 



,{p) = S{p~ Po) . 



In this case the particle distribution stays the same for all 
times and the expectation values E{i?mcomp) — po and 
E(T) — E{Tinit) = N/n are constant. The appropriate 
translation between the mass density PDF and the par- 
ticle distribution is a convolution which maps the S-peak 
to the Poisson-distribution and matches the expected val- 



ues: 



f(fc) = {R®Px){k) 



= P. 



= P.nTlSk)R{p)Ap. 



(2) 



It is easy to see that this convolution works for the in- 
compressible case. From that it is suggesting to take the 
same convolution for compressible turbulence with the 
measured R[p) instead of the delta peak in the incom- 
pressible case. In this way the convolution generates a 
discrete particle distribution from the continuous distri- 
bution of density. In Fig. [3] we show the PDFs for run 
a and c (The PDFs for run b agree similarly well, but 
are omitted here). The solid lines T{k) are the result of 
the convolution of the measured R{p) with the Poisson 
distribution. It reproduces in all cases nearly exactly the 
measured particle distribution T{k). 

In order to obtain an analytical expression, we use the 
relation between the deviation a of p and the r.m.s. Mach 
number M for isothermal supersonic turbulence intro- 
duced by Passot and Vazquez-Semadeni [ij] and Nord- 
lund and Padoan [l^, [l^ : 



This translates to the relation 



£2 (R) 



(3) 



of the standard deviation as of In p and the Mach number 
M. 

If we convolve this lognormal distribution with a 
Poisson-distribution using eqn. ([2|) we obtain 



f{k) = 

with 

E = 



1 



27rcrs fc! 



p^ ^exp 



2ai 



(4) 



E(r) 
MR) 



Hence, we can describe the particle-distribution with only 
one fit parameter /?. In Fig. [3] we show the agreement be- 
tween the particle-distribution and this theoretical PDF. 
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FIG. 3: The convolution operator maps the mass density PDF 
R{p) on the particle distribution T{k) in good agreement with 
the measurement T{k). Top: run a, bottom: run c 



This figure was obtained with a value of ~ 0.37. This 
value lies somewhere in between the value /? — 0.26 of 
the recent simulations of Kritsuk et al. \1TV a nd (3 — 0.5 
of the MHD simulations of Padoan et al. [l6| . 
Structuring time After this findings we can determine 
how long it takes for the particles to form global struc- 
tures. This time tg is equivalent to the time in which 
the particle PDF goes from the initial Poissonian to the 
final distribution given in eqn. (U). To estimate this time 
we define an operator A[r] on the particle distribution T 
which assigns the initial Poisson distribution T = P{k) 
to 1.0 and the final stationary distribution T = T{k) to 
0: 



- Ep) dp ((1 - • ^ (^•) + ^[^] ■P{k)-T [k) 



Min 



where the value A[T] is obtained from the minimization 
procedure. A plot of the temporal development of A[r] 
is shown is Fig. [51 

To give an estimate for the time scale ts of the particle 
clustering, an exponentially damped oscillating function 
is fitted to the curve A[T](t) with decay proportional 
to exp(— i/ig). In Table HIl values for this structuring 
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FIG. 4: Particle distribution averaged over several time steps 
T (k) and the theoretical distribution T (k) obtained horn 
equation Q with 13 = 0.37. 
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FIG. 5: Transition of the particle distribution from the Pois- 
sonian to the T'(fc)-distribution characterized by A[r](t) and 
compared to the exponential fit f{t). 

time ts as well as the dynamical time t^yn are given. 
The comparison shows that the structuring time is sig- 
nificantly shorter than the dynamical time, although the 
ratio ts/tdyn grows slightly with increasing Mach num- 
ber. 



TABLE II: Comparison between the dynamical time tdyn and 
the structuring time ts 



Run 


Mach 


^dyn 


ts 


a 


0.4 


7.8 


1.0 


b 


1.4 


2.2 


0.33 


c 


4.6 


0.62 


0.11 



Conclusions and outlook: In this Letter we have estab- 



lished a direct relationship between the particle distribu- 
tion and the PDF of the density fluctuations in isother- 
mal compressible turbulence. The analytical expression 
for the particle distribution confirms the relation between 
the standard deviation of In p and the mean Mach num- 
ber. The proportionality constant /3 was obtained from 
the particle distribution. The Lagrangian viewpoint of- 
fers further studies on the formation of density struc- 
tures. In this paper, we studied only the low order statis- 
tics. The next natural step is to consider higher order 
moments of the particle number in a given ball of ra- 
dius r. A first step in this direction for an unrealistic 
compressible Kraichnan flow was done in Bee et al. [l^ . 
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